home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / CUJ9204.ARJ / 1004030A < prev    next >
Text File  |  1992-06-02  |  2KB  |  70 lines
  1. /* Listing 3. */
  2. typedef struct {
  3.     float cos1, sin1, cos2, sin2;
  4. }      ARG_FF_2;        /* vector 4 pairs */
  5. ARG_FF_2 cosf_sinf_2(x1, x2)
  6.     union {
  7.     float flt;
  8.     int iflt;
  9.     }     x1, x2;
  10. {                /* 2 pair single precision
  11.                    sin/cos function */
  12. #include "float.h"
  13. #if FLT_ROUNDS != 1
  14.     #error "rounding mode not nearest, fix code"
  15. #endif
  16. #include <errno.h>
  17. #ifndef errno
  18.     extern int errno;
  19. #endif
  20. #include <math.h>
  21. #define M_2_PI        0.63661977236758134308
  22. #define T2PI        2*M_2_PI
  23. #define TP2        T2PI*T2PI
  24. #define TP3        TP2*T2PI
  25. #define TP4        TP2*TP2
  26. #define TP5        TP3*TP2
  27.     ARG_FF_2 res;
  28.     double tn, td, r;
  29. /* integer compare with 0 same as float, for IEEE
  30. ** since arg comes in int register, this is faster
  31. **
  32. ** doing everything in double, we won't lose accuracy
  33. ** by converting arg to multiple of PI/2
  34. ** this allows range reduction by subtracting an integer
  35. **
  36. ** reduce to range +- 2, divide by 2 later
  37. ** when it cannot underflow */
  38.     tn = x1.flt * M_2_PI + (td = x1.iflt >= 0 ?
  39.             4 / LDBL_EPSILON : -4 / LDBL_EPSILON);
  40.     if (fabs(x1.flt * M_2_PI) >= 4 / LDBL_EPSILON |
  41.     fabs(x2.flt * M_2_PI) >= 4 / LDBL_EPSILON)
  42.     errno = ERANGE;
  43.     tn -= td;
  44.     tn = x1.flt * M_2_PI - tn;
  45.     td = tn * tn;
  46. /* divide arg by 2 and rationalize numerator and denominator
  47. ** numerator of rational approx for tan(x1/2)
  48. ** Horner polynomials 1st time */
  49.     tn *= 886.77348 * TP4 + td * (-99.398954 * TP2 + td);
  50. /* denominator */
  51.     td = 886.77346 * TP5 + td *
  52.     (-394.98971 * TP3 + td * 14.425694 * T2PI);
  53. /* cos, sin half angle formulae, rationalized */
  54.     res.cos1 = (td * td - tn * tn) / (td * td + tn * tn);
  55.     res.sin1 = (tn * td + tn * td) / (td * td + tn * tn);
  56. /* copy 2 */
  57.     tn = x2.flt * M_2_PI + (td = x2.iflt >= 0 ?
  58.             4 / LDBL_EPSILON : -4 / LDBL_EPSILON);
  59.     tn -= td;
  60.     tn = x2.flt * M_2_PI - tn;
  61.     td = tn * tn;
  62. /* distribute terms to finish polynomials quicker */
  63.     tn *= 886.77348 * TP4 - td * 99.398954 * TP2 + td * td;
  64.     r = 886.77346 * TP5 - td * 394.98971 * TP3;
  65.     td = r + td * td * 14.425694 * T2PI;
  66.     res.cos2 = (td * td - tn * tn) / (td * td + tn * tn);
  67.     res.sin2 = (tn * td + tn * td) / (td * td + tn * tn);
  68.     return res;
  69. }
  70.